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Abstract. A self-consistent calculation scheme for correlated electron systems is 
created based on the density- functional theory (DPT). Our scheme is a multi-reference 
DPT (MR-DPT) calculation in which the electron charge density is reproduced by 
an auxiliary interacting Pcrmion system. A short-range Hubbard-type interaction is 
introduced by a rigorous manner with a residual term for the exchange-correlation 
energy. The Hubbard term is determined uniquely by referencing the density 
fluctuation at a selected localized orbital. This strategy to obtain an extension of the 
Kohn-Sham scheme provides a self-consistent electronic structure calculation for the 
materials design. Introducing an approximation for the residual exchange-correlation 
energy functional, we have the LDA+U energy functional. Practical self-consistent 
calculations are exemplified by simulations of Hydrogen systems, i.e. a molecule and a 
periodic one-dimensional array, which is a proof of existence of the interaction strength 
[/ as a continuous function of the local fluctuation and structural parameters of the 
system. 



I This paper is to be published in J. Phys. Condens. Matter (2007). 
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1. Introduction 

Inclusion of tlie sliort-range correlation effect (SRCE) is a long-term request for the 
first-principles electronic structure calculation based on the density functional theory 
(DFT).[Tl [2] In principles, it is possible, since the strategy introduced by Hohenberg, 
Kohn and Sham was shown to be given by a rigorous variational principle. [3], HI El E] 
Although the method should give formally an exact calculation scheme for the Coulomb 
system, the energy density functional is not perfectly known at present. Plausible 
approximation schemes have been proposed and utilized. [21 [71 El [9] However, they have 
their own limitation. For example, the local-density approximation (LDA) is known to 
conclude a metallic ground state for the Mott insulator La2CuO4.[T0l [TTl [121 [131 [IH 
[T5l [T6] This failure of LDA is a central problem of DFT for which we hope inclusion of 
SRCE to be a solution. Especially, when LDA gives near degeneracy in the ground state, 
proper treatment of SRCE can lift the degeneracy to have the non-degenerated ground 
state implying formation of the Mott gap. This assumption may be widely accepted as 
a natural conclusion according to the study of the Hubbard models. [TTl [18] 

Here we should note that the Kohn-Sham scheme has flexibility and it can be 
adjusted even for the Mott insulator. If we introduce a wavefunction of an entangled 
state as the Kohn-Sham ground-state wavefunction, the excitation spectrum for the 
Kohn-Sham system may change. This implies that response of the system has changed. 
Considering the adiabatic shift of the ground state as a function of some outer 
parameters like the external electro-magnetic field, there should be an essential change 
as a consequence of the introduction of SRCE in the Kohn-Sham scheme. Even if we 
consider the density functional theory for the ground state of the Coulomb system, 
this extended scheme allows us to handle the correlated electron system by the density 
functional theory. 

Thus we have yet many possible approaches for the practical computation as 
realization of the Kohn-Sham scheme in an extended formulation. Actually, the Kohn- 
Sham equation is regarded as an auxiliary equation to realize the optimization process 
of the single particle density n(r). In this paper, we consider this physical quantity as 
a central order parameter of the electron system. Usually, a system of non-interacting 
Fermions is utilized to describe n[r) in the Kohn-Sham scheme. Interestingly, we are 
allowed to consider interacting Fermion systems, which can be used to replace the non- 
interacting Kohn-Sham system. This is called the multi-reference density functional 
theory (MR-DFT).[T9l [201 [22l [231 [25] To develop direct description of a Mott insulating 
state, one of the authors defined a kind of MR-DFT.[26] Utilizing this formulation 
called the extended Kohn-Sham scheme (EKSS), one has a chance to detect Coulomb 
suppression of fluctuation, which is not found in n(r). 

The interacting Kohn-Sham system has been originally motivated in the hybrid 
approach with the configuration interaction (CI) scheme in the theory of the quantum 
chemistry, [m [201 [211 [22l [23l [241 [25] In the hybrid density functional theory, people 
utilized 1) a full or a part of elements or integrals of the density matrix [22] or 2) 
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restriction of the searching space [23] in the constrained minimization to define the energy 
density functional Knowledge on the modified energy density functional, however, are 
not enough. A proof of existence of the minimum in the constrained search is demanded. 
On the contrary, it is possible to formulate MR-DFT in another way by referring the 
original Levy-Lieb energy functional. [261 EZj 

In this paper, focusing on the fiuctuation reference method, [27] we will discuss 
a self-consistent calculation scheme of MR-DFT. The method is shown to be a kind 
of the renormalization method to find a fixed effective interacting Hamiltonian. A 
practical approximation for the residual exchange-correlation energy functional allows 
us to confirm that the scheme do give the self-consistent solution. We will give a 
report on the first application of our scheme in two types of the Hydrogen systems. 
If we introduce a local density approximation after replacing the residual exchange- 
correlation energy functional by the ordinal exchange-correlation energy functional, the 
obtained energy functional is a kind of the LDA-I-U energy functional. However, our 
approach is different from the former LDA-I-U approaches, [281 EH], [30] because we follow 
the fluctuation reference method and not primarily looking at the excitation spectrum. 
Clear difference from the LDA-I-U approach can be seen in the fact that we are able 
to avoid the clued approximation replacing the residual exchange-correlation energy 
functional by the ordinal exchange-correlation energy functional. 

The structure of the paper is as follows. In Sec. 2, we introduce our energy 
functional. The functional is a wave-function functional. The variational principle is 
shown. In Sec. 3, the idea of the fluctuation reference is introduced. The uniqueness 
theorem of the U term is briefly reviewed. We discuss the extended Kohn-Sham 
Hamiltonian as a flxed point Hamiltonian in MR-DFT in Sec. 4. In Sec. 5, importance 
of the density fluctuation to determine the correlated nature of electron systems is 
discussed. In Sec. 6, we introduce a practical application of the method with two 
Hydrogen systems. Final discussion and summary is given in Sec. 7. 

2. Energy functional 

We review the formal theory of the extended Kohn-Sham scheme (EKSS).[26,j We 
consider a non-relativistic electron system with electrons in an external scalar 
potential fext(r). The Hamiltonian operator that we consider is. 




(1) 



The kinetic-energy operator is given by. 




and the inter-electron repulsion is. 



V;^ = 1 / £r£r' 
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The ground state |\1/gs) of the system exists and gives the lowest energy Eq and the 
single particle density as, 

^0= (^GslHcl^Gs). (2) 

nGs(r) = (^Gslri(r)I^Gs)- (3) 
Here n(r) = Xlo- V'(l(^)V'o-(^) with the electron field operator ipai^^) satisfying 

We know the following density functional theory. |5j For a normalizable wavefunction 
\E' with a finite kinetic energy, the single particle density n(r) of \I' and |V(n(r))^/^p are 

in a set of integrable functions in M^. A set is a set of functions / for which J p 
and y |V/|^ are finite. We consider a minimization scheme with respect to n(j) > 

such that n(r)^/^ G H^i^') and ^ n(r)d^r = N. This class of functions is called In- 

Since a minimizing sequence of a positive quadratic form in Hi{M.^'^) has a limit, 
and since the Harriman construction |3 11 15] ensures existence of \E'' giving n(r) G Tn, one 
can introduce a universal energy functional F[n\ which is called the Levy-Lieb energy 
functional and defined by 

F\n\= min (^'|f + Kel^') • (4) 

<I"^n(r) 

Utilizing this energy functional, we can construct the minimization process of EKSS. 
To formulate it, let us consider a set of orthogonalized normalizable functions {(f)i{r)}, 
the creation and annihilation operator c]^ and Cjo-, and a number operator fiia- = c\^Cia- 
with respect to 0i(r). Expectation values njo- = {^\nifj\^) are given for a state |\1'). We 
introduce another density functional. 

Fu[n]= mill ($'|r+^ V(n,|+ni^-%-ni|)2|^f'). (5) 

i 

There is a minimizing state for any nij) G X/y. 

As the ordinal Kohn-Sham scheme, EKSS ensures that the total energy E^ and 
the single-particle density nGsli") of the ground state are reproduced. This is due to 
the definition of the optimization process utilizing the Levy-Lieb energy functional. The 
physical phase space of |\E') is divided into pieces specified by their single particle density 
n{j). Then, the minimization process is decomposed into the inner process with respect 
to within the subspace given by n{j) and the outer process with respect to n(r). 

If we further make an attention on the Hadjisavvas-Theophilou scheme, [6] we can 
show EKSS in a rigorous manner. This process is easily shown in the next equality. 



= (^Gs|T + Kel^Gs) + / nGs(r)i;ex.t(r)rf= 



min<^ min (v^jT + Kel^E') + / n(r)t;e^j(r)(iV 

*-^n(r) 
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= mill I mill (\I>'|T + ^ "^(^it + ^ii " " 

n I *'^n(r) Z ^—^ 

- Fu\n\ + ^ n{j)v^^t{j)d^r\^ 



mm < mm 



2 / r — r' 



lin J (xl/'lf + f/ ^ Ti^rn^ + I 5^ 



= miiiGc/l^']. (6) 
Here, n^, is the density associated with 
n^{j) = (^|n(r)|^) . 

Thus we have found that the minimization process of a wave-function functional (5;7[\[^'] 
gives the exact value of the total energy of the system. 

In a general form, the energy functional Geks[^] of EKSS is given in the next 
formula. 

GEKsm = + Kcdl^) - min (^'|f + Kedl^') 



+ F[n^] + J d^rvc^t{r)n^{r 

2 J |r — r I 



+ E,^c[n^] + j c/Vfext(r)n,i,(r) . (7) 

Here the operator V^ed denotes a generalized operator counting fluctuation or hidden 
order parameters which are in a higher order than that of n(r). The operator 
has to be a positive semi-definite and be bounded from above. When minimizing 
Geks\^] with respect to which is an auxiliary wavefunction, the value of Geks\^] 
becomes £"0. This is easily seen by looking at the first line of Eq. ([7j), in which 
(\I/|T -|- V^edl^) — min^/^n^ (\E''|T + V^redl^') > becomes zero at the minimum. At this 



A first-principles calculation for correlated electron systems 



6 



minimum point, \1/ gives the minimum value of the expectation value {"^IT + Kedl^) 
within a phase space of wavefunctions whose single particle density is n^,. Now, the 
density functional -F[n^] + J (PrVextiT')'n^i^) becomes minimum, when n^(r) is equal to 
the single-particle density of the true ground state nGs(r). Thus, the total minimization 
is achieved, only if n,i,(r) = nGs(r) and if gives the minimum of (\E'|T-|- V^edl^) within 
the phase space of wavefunctions which give nQs{r). 

One would find that Eq. ([7]) is nothing but the definition of i?rxc[^*]- Formally, V^ed 
is arbitrary, since redefinition of Ej-xdn^] keeps the equality. Moreover, the kinetic term 
and the Hartree term are not necessarily given by the formula in Eq. ([7j). At present, 
we just follow the conventional idea that the Hartree-type approximation would close to 
the answer, when we know a priori the density n{r). Using the usual Kinetic energy of 
Fermions with the electron mass, we have determined Ej-^dnq,]. This guideline may be 
explained in the following manner. If we know that n(r) is the proper order parameter, 
it would be natural to expect that the explicit energy functional written in n(r) with 
the Hartree term reflects dependence on the structure of the materials at the first stage. 
The electron charge density acts as a source and creates the scalar Coulombic field. In 
addition, introduction of the Fermion kinetic energy (\E'|T|\E') keeps the system from 
the collapse to the Bosonic solution. The reason why we conclude the above statement 
is that the variable of the theory is n(r). The Kinetic energy functional, however, has 
another meaning as discussed in Section [71 

An important point for the density functional theory is that we can find continual 
refinement for the improvement. Introduction of (\l/|V^ed|^) shifts the energy functional 
so that |\E') represents a correlated electron state. Using the entangled state, expression 
of the energy functional is modified. In the new description, explicit evaluation of the 
energy is done with the Hartree term, the kinetic energy and the fiuctuation. If the 
residual correlation energy functional i?j.xc['^i'] becomes small in its ratio to the total 
energy by this modification, we notice that the fiuctuation has emerged. Now we start 
to explain the idea in detail. 

To proceed, we need to consider functional differentiability. [9] For this purpose, all 
of the energy functional defined above should be replaced by the Legendre transforms 
of them. The technique was introduced by Lieb.[5] To specify the problem, we consider 
(5[7[\E']. By making a variation with respect to (^|, we have an extended-Kohn-Sham 
equations (EKSE). 

i 

+ J2^il-2n,)Y,n.,.m = m)- (8) 

i a 

Here fii = X^o-^^.o"- ^ Lagrange multiprier E is introduced to keep the norm of |\I') to 
be one. Here the effective single particle potential fefr(r) is given by. 




T + Vcs{r)n{r)d^r 
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The charge density n(r) is given by 

niv) = j2mii^)Mrm- (10) 

Please note that we have not yet given a determination metliod of {(f)i{r)}, but that the 
variational principle holds always rigorously. 

We solve the auxiliary one-body problem given by as, 

-— Ar + t;eff(r)|xKr) =£iXKr), (11) 

in which xK^') determined to be normalized and orthonormal. If we construct a 
set of creation and annihilation operators d]^, di^a associated with Xi{^)-, the effective 
many-body problem is found. 

K l,a i i a J 

Note again that fii^^ = cl^Ci^a is defined by (/)i(r). In a crystal, the index I may be a 
combined index of the crystal momentum k and the band index n. One may call EKSE 
defined by Eqs. ( ITTl) and ( |T2i) a first-principles Anderson model or a first-principles 
Hubbard model. 



3. A comment on the uniqueness of the model 

In principle, EKSS works irrespective of the form of V^ed, if we can check existence of 
the minimum of (\E'|T + K-ed|^) and its bound. This fact tells us about flexibility of 
the theory. A big class of effective Hamiltonians exists and each auxiliary system is an 
extended Kohn-Sham model. Thus, we need to have a rule to select a properly chosen 
effective model for a practical calculation. In other words, there should be a guiding 
principle to determine (j(/[^E']. The idea is that there has to be a physical quantity which 
is in a higher order than ?T,(r) and specifies the model. 

At the beginning, we need to understand nature of G;7[\E'] to construct the best 
fitted model. To make the discussion concrete, let us consider a U term in our theory. 
For a given normalizable localized orbital 0i(r), density fluctuation is determined as 
follows. 

(li) = (("-i.T + - ^i,T - ^i,^?) . (13) 

A key observation is that the fluctuation counted by the model may be suppressed, if 
the minimizing ^ changes when the value of U is increased in eq. (fT2l) . 
The U term in G[7[\I''] is given by the next energy functional. 

(^iKedl^) = ^{^llim , (14) 

A requirement is that the U term has to be bounded from below and from above. This 
is guaranteed in the above expression, since the quadratic form is positive-semi definite 
and the lemma below holds. [27] 
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Lemma 1 

(n^) is real. The next inequality holds. 

< {ri) < 1 . (15) 

We also have next few statements, which are given without proof here. 

Lemma 2 Assume that the ground state of a Coulomb system given by fext(r) exists. 
i) The ground state of a corresponding extended- Kohn- Sham model Gu\^] with 
a given positive U exists, ii) For fixed n{r), F{U) = min,i,_,„(\E'|T + ^ri?|\l/) is a 
continuous function of U. iii) // a state is the ground state of Gui[^] and Gu2[^] 
with < Ui < U2 simultaneously, is the ground state of Gu[^] in a finite range 
[UuU2] ofU. 

The proofs are given in another paper. [27] Finiteness of (n^) is utihzed for the proof of 
the continuity. The constraint for the degeneracy of the Coulomb system is not required 
in Lemma [2j 

If we increase U from zero, the effective interaction in Eq. ([H]) brings the system in 
a correlated regime. The change results in the suppression of the fluctuation considered. 
Thus, the U term can control the value of {nf). For the original Coulomb system, we can 
also determine (11^)03 in principle, once we fix 0i(r). We are thus allowed to compare 
the fluctuation of the original system and the extended Kohn-Sham system. There could 
be an adjusted value of U for which (n^) of EKSS is identical to that of the Coulomb 
system. 

At a first glance, this point is not so important, since the density-functional theory 
tells nothing about fluctuation or correlation functions. The Kohn-Sham wavefunction 
is introduced to determine the minimization process with respect to n(r) and do not 
have direct relevance in itself. However, if the given extended Kohn-Sham system is 
properly written in a multi-reference description, and if the obtained extended Kohn- 
Sham model reproduce an essential nature of the original system, the theory may have 
gone beyond the original concept of the density functional theory. 

For example, introduction of U can make the extended Kohn-Sham system the Mott 
insulator. The solidification caused by suppression of the density fluctuation given by 
(zZii) niay be detected in practical calculation. As discussed in Sections E] and [3, we can 
judge whether the system is the Mott insulator or not. Thus reproduction of important 
fluctuation can be a key procedure to have a good description of some materials. 

In a previous work, Kusakabe has shown a statement on uniqueness of the U term. 
We have the next exact statement. 

Theorem 3 Assume that the ground state of a Coulomb system is non- degenerate. A 
proper extended- Kohn- Sham model given by Gu\^] which has a non- degenerate ground 
state and reproduces both riGs(r) and (nf)GS uniquely determined, or it does not exist. 

This is a principle of our fluctuation reference method. 

The restriction on the degeneracy of the Coulomb ground state is strict in the above 
theorem. Some systems are known to have degeneracy in the ground state. As for the 
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degeneracy due to the spatial symmetry, the condition may not be a problem, since 
we are allowed to consider an outer scalar field which breaks the symmetry. Internal 
symmetry considered in the present description of the many-electron system with Eq. 
([1]) is the electron spin. We may have degeneracy due to the internal symmetry of this 
spin degrees of freedom. As for the trivial degeneracy coming from the SU(2) symmetry 
of the total spin, an external magnetic field will lift the degeneracy via the Zeeman 
splitting. If we change the structure of atomic configuration, effective spin interactions 
in the system change to lift the degeneracy in some cases. 

4. Renormalization of the extended Kohn-Sham model 

We now clarify that the self-consistent determination of the extended Kohn-Sham model 
is a sort of the renormalization process. We consider Eq. ([8]) or Eq. ( fT2l) . The set of the 
solutions of Eq. ( ITTi) can be used to create <pi{r). In each step in the self-consistent 
loop, n{r) is changing gradually and thus too. What can be fixed in the process 

is an algorithm to make 0j(r) from Xii''^)- 

More precisely, considering a lattice structure, we can diagonalize the single-particle 
part by Bloch waves xK"^) = Xn,k(i")- The orbital is at first specified by a combined index 
/ with the band index n and the crystal momentum k. A unitary transformation from 
the Bloch states to the Wannier states may be useful to define 0i(r) as 0m(r — R-j)- We 
suppose that i denotes an m-th localized orbital at a Wannier center Rj.[32l |33] If we 
fix the selection of the relevant bands to create the Wannier states, the self-consistency 
loop to find a solution of Eqs. ( ITTi) and ( |T2i) is well defined and it may converge. 

In the model of Eq. ([8]), the scattering channels given by the effective interaction 
term are open only within a subset of xK^)? which is determined by the selection of 0i(r). 
In other words, is expanded in dj^ in a specified n-th band only. The scattering by the 
U term is restricted within this band and no direct interaction with other bands exists. 
Thus the definition gives a separable form of the effective interaction. If scattering 
processes due to the effective interaction are completely restricted within selected bands, 
the form is called separable. 

If the effective interaction is written in terms of the field operators V'o-lr), and if 
the interaction strength g{r,r') is not written in the separable form, there should be a 
finite amplitude for the scattering channel from one band to all the other bands. Thus, 
to solve obtained EKSE is as hard as the original Coulomb problem. But if the Fermion 
scattering processes due to the effective interaction are restricted in a specified sub 
space of the whole phase space, reduction in the many-body description is achieved. If 
relevant scattering processes are properly chosen in the effective model, and if the total 
self-consistency is achieved, the obtained Hamiltonian should be a fixed Hamiltonian. 
The point is that the orbitals to describe the effective interaction have to be determined 
self-consistently. 

Arbitraryness of 0j(r) actually allows us to have the fixed Hamiltonian. We can 
redefine the U term in an optimization process of G[/[\E'] by making use of 0i(r) given by 
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the selected n-th band in the calculation. If 0i(r) given in a step of the self-consistency 
loop is not perfectly expanded in the former set of wavefunctions in the n-th band of Eq. 
( |TT|) . we can reconstruct 0i(r) as a new Wannier orbital in the next step starting from 
the obtained n-th band. This approach to redefine the effective interaction is regarded 
as a renormalization process. The final fixed-point Hamiltonian would be described 
in a specified relevant sub-space whose dimension is much smaller than the original 
problem. Irrelevant scattering processes are smeared out from the theory. As for the 
electronic charge density n(r), which is an essential quantity to determine the structure 
or the atomic configuration of a material, introduction of the renormalization process do 
nothing harmful, since the obtained effective Hamiltonian gives the ground state charge 
density and the ground state energy. 

5. Density fluctuation 

Density fluctuation (n?) plays an important role in our theory. The reason why we select 
this quantity as a physical quantity second to n(r) may be explained as follows. 

This quantity has a value depending on the environment around 0i(r). Consider a 
d orbital of a cupper atom as an example. The fluctuation on the orbital would be not 
small, when cupper atoms form a bulk metal. But, if the atom is in a cupper oxide, the 
fluctuation should be reduced on it due to SRCE. 

In an ideal case, we can have a correlated electron state as the ground state, whose 
electron density n(r) is the same as another uncorrelated state but it has a different 
fluctuation on 0i(r). The theory in Section [3] tells us that an effective many-body system 
properly describing both n(r) and the fluctuation (n^) on 0i(r) is uniquely determined, 
if it exists. The ground state of the model would have a correlated state and sometimes 
it becomes even the Mott insulator. A typical example may be the Heitler-London state, 
which is an entangled singlet state. 

Considering both the uncorrelated metallic state and the entangled state in a 
correlated regime, we can easily understand the essential behavior of (n?) as follows. 
For a nearly uncorrelated metal, it is easy to show that (n^) = 0.5. However, it should 
be zero for the Heitler-London state of the Hydrogen molecule, as exemplifled in Section 
El 

We may deflne the Fermi level Ep for convenience, once Eq. ( |TT1) is solved with 
a flxed number of electrons. The wavefunctions X;(r) is grouped in bands. For each 
band, a unitary transformation to a localized orbital 0i(r) is given. The typical value 
of fluctuation on it is classified in the next lists. 

(i) If (pi{r) is deep below Ep, (n^) = 0. This is because the orbital is doubly occupied. 

(ii) If 0j(r) is far above Ep, {nj) = 0. This is because the orbital is empty, 
(in) If 0i(r) is around Ep and if the state is uncorrelated, (n^) = 0.5. 

(iv) If 0i(r) is around Ep and if the state is correlated, (n^) = 0. 
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We have to select 0i(r) to keep symmetry of the system, otherwise we will encounter 
difficulty in description of the system. Another important comment is that, if we choose 
an extended wavefunction as 0i(r), the fluctuation on it may approach to (n^) = 1 in a 
correlated regime. 

6. Determination of U in the Hydrogen Systems 

In this paper, we consider Hydrogen systems to demonstrate that it is possible to 
determine 1) the self-consistent solution of the extended Kohn-Sham scheme, and 2) 
the interaction parameter U, in practical calculations. Since the relevant orbitals are 
only Is orbitals in the Hydrogen systems, the electronic structure is easily tractable. We 
select two systems, i.e. the Hydrogen molecule and a one-dimensional lattice structure. 
(Figure 1) The former example shows that an entangled state is obtained as a self- 
consistent solution of the extended Kohn-Sham model. The U term is determined by 
fitting the local fluctuation of an accurate CI calculation for the Hydrogen molecule. The 
latter seemingly artificial configuration of a Hydrogen chain with a periodic boundary 
condition is introduced to show that a Mott-insulating state is obtained as a self- 
consistent solution. 

For both of these systems, the extended Kohn-Sham equation is given in Eq. ([8]). 
The value of U is identical for every site indexed by i, because of the symmetry of the 
system. More precisely there are the C2 symmetry (the mirror symmetry with respect 
to the center of the molecule) for H2 and the translational symmetry (invariance for 
uniform shift by the lattice constant a) for the chain. For both of the system, we have 
no spontaneous symmetry breaking causing the charge density wave, because the final 
solution of EKSE is non-degenerate. 

6.1. Self- consistent calculation method 

The self-consistent calculation is realized by adopting an algorithm with two nested 
loops. The outer loop is the determination of the CI configuration of the effective 
many-body problem. The inner loop is the diagonalization of Eq. (fTT!) to obtain Xii''^)- 
The index / is / = 1, 2 for a bonding state and an anti-bonding state in the molecule. 
While it is / = (n, k) with n = 1 and k = 0, ■ ■ ■ , N — 1 for the chain, n = 1 corresponds 
to the Is band. We define 0j by 



01 = (Xi + X2) , 

02 = -7^ (Xl - X2) , 



for the molecule and the Wannier state 
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Figure 1. The calculation cell of the Hydrogen systems, (a) the Hydrogen molecule 
and (b) a Hydrogen chain. The inter-atomic distance is R[A] or a[A\ for the molecule or 
the chain. The system in (b) consists of 10 atoms with a periodic boundary condition. 
The outer cell is for the many-body calculation. The inner cell denoted by dashed lines 
is for the determination of the single-particle orbital xi (r) . 



for a chain with atoms. Note that the size of the outer cell in the x direction is Na. 
Xi,k is the Bloch wave in the first Is band with the crystal momentum p = 27Tk/{Na) 
in the chain direction. Xi = ai is the x-coordinate of the i-th atom. (Figure 1) In the 
present systems, we can determine the transfer matrix element by 

tij = J 0:(r) |-|^Ar + i;eff(r)| 0,(r)dr. (16) 

Here, a dependence does not appear because the system is non-magnetic. We select a 
typical transfer energy to as that between the nearest neighbor pair of orbitals. The 
U term is then introduced and Eq. f[T^ is diagonalized. For the case of the chain, we 
utilize the numerical diagonalization with the Lanczos algorithm. Here, the problem is 
solved for a fixed U/to. Fixing the CI configuration, the one-body problem of Eq. (fTTll 
is solved self-consistently. Then, using the determined new Xh the effective Hubbard 
model is rebuilt. The self-consistency on the CI configuration is checked in the outer 
loop. Actually, we can reach the totally self-consistent solution. 
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The residual exchange-correlation energy functional is rewritten as follows. 

E,^^[n^] = F[n^] - min {^'\f + V^^') - \ f "^^P^^Srdh' 
^'^n^ 2 J |r - r'l 

= FM - min (<i.'|r|<i.') - - [ I^iMl^dVdV 

^ ^ *'^n*^ ' ' ' 2 J \r-r'\ 
+ min (<l>'|f 1$') - min (*'|f + Kodl^') 

= ^xcM + min ($'|f 1$') - min ($'|f + Kcdl^')- (17) 
One way to treat the above expression is utilizing the next approximation. [51] 
min (^'|f + l/redl^') ^ min (^'|f 1^') + min — 

i 

= min (^'ITI^'). (18) 

If the serching space of $' in Eq. (fT7|) is the set of the single Slater determinant 0', and 
if min<^/_,„^(0'|T|0') = min^,/_^„^ (\E''|T|\E''), E^dn^] is the same as the ordinal exchange- 
correlation energy functional. Note that and $' are multi-Slater determinants. This 
is true if we consider the Legendre transform of each expression. If Eq. (|T8|) is adopted, 
we see that i^rxcf'^*] — -Exc[^*]- Then, i^^xc[^*] is approximated by the local-density 
approximation. [36] The treatment of Eq. (fTTll will be reconsidered in Section [71 For 
the actual calculation in the inner loop, we utilized the plane-wave expansion technique 
with the soft pseudo potential. [37] To use the pseudo potential with LDA does not harm 
the purpose of the present MR-DFT calculation, which is planned to show existence of 
self-consistent solutions. The cut-off energy is set to be 40 [Ry]. The conjugate-gradient 
technique is used to optimize the Kohn-Sham orbitals X/(i")- The actual calculation 
was done using a computation code called ESopt, which was originally developed by T. 
Ogitsu and maintained by K.K. 

6.2. Reference calculation 

As the reference calculation, we refer to the result obtained by the complete-active-space 
configuration-interaction (CASCI) theory [T9] for the Hydrogen molecule. The Gaussian 
basis set is utilized. The CAS wavefunction is prepared to incorporate all the resonating 
features arising in the H2 molecule. Another MR-DFT approach, the CASCI density 
functional theory (CASCI-DFT) was also examined. In the CASCI-DFT calculation, 
the CI configuration is taken from the CASCI calculation. The detailed description on 
the exchange-correlation energy functional used in CASCI-DFT is seen in Ref. |25j . 
The fluctuation on the Is orbital is obtained as a function of the inter-atomic distance 
as shown in Figure 2. 

When the inter-atomic distance is the equilibrium value R = 0.740A, the fluctuation 
is close to 0.5. This result tells us that the two-electron system of the Hydrogen molecule 
is in a weak correlated regime, when the system is in equilibrium. However, when R 
becomes larger than lA, the fluctuation is rapidly suppressed. This is seemingly natural. 
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Figure 2. The density fluctuation (nf ) at a Is orbital of the Hydrogen molecule. For 
the inter-atomic distance (n?) is obtained by the CAS-CI calculation (circles) 

and the CAS-CI-DFT calculation (crosses) . 



since the system should approach the Heitler-London hmit when R ^ 0.740A. Crossover 
region is thus shown to be i? ~ 2A. 



6. 3. EKSS calculation of the Hydrogen molecule 

The MR-DFT using the extended Kohn-Sham scheme is apphed to the Hydrogen 
molecule. [35] Formally, the value of the fluctuation should be given as a function of 
i?[A] and t/[Ry]. However, we obtained U for given R[a\ and (n^). (Figure 3) In this 
case, fixing (n?) is equivalent to fix f/ = f//to- The value of to is given, when the inner 
loop is converged. The value of f/ = [/to is thus known after finding a self-consistent 
solution. The solution is obtained for each fixed U and R. By comparing (n^) of the 
effective model with that obtained by CASCI, U is determined uniquely. (The solid line 
in Figure 3) 

Since we utilize the pseudo-potential method, the obtained 0i(r) in the model is 
not the same as that given by CASCI. Thus the estimated value is an approximated 
one. In principle, evaluation of (n^) using 0j(r) in CASCI is possible. An important 
point is that the obtained value of U changes continuously and monotonously. Thus, in 
this numerical evaluation, determination of U is possible. 
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Figure 3. The effective interaction parameter U obtained by the extended Kohn- 
Sham calculation for the Hydrogen molecule. The dashed lines are the values of U 
for R = 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.8, 2.0, 3.0[A] from the top to the bottom. The 
density fluctuation (n?) is counted on a Is orbital of a Hydrogen atom. By adjusting 
iHi) to the result in Figure 2 for each value of R, optimized U is obtained. (The solid 
line) 

6.4- A one- dimensional Hydrogen array 

As the second test calculation, we consider an array of Hydrogen atoms. The 
configuration is imaginative, since the structure is not stable and inter-atomic forces 
are finite. But, to consider a simple Mott insulator, this artificial configuration is very 
useful. 

We consider a periodic boundary condition with 10 atoms {N = 10) in an outer 
simulation cell. (Figure 1) Since the system does not show spontaneous charge ordering, 
electron charge density n{r) has the periodicity that is same as that of the array. Thus, 
we can consider an inner unit cell containing a single atom in it. Within the second unit 
cell, 77, (r) is kept in the simulation. 

For a multi-reference state, we have an expansion. 

= J2C^\^a), (19) 

a 

1^")= rK^tiidio)- (20) 

m=l n=l 

Here, a is an index specifying the CI configuration. Considering up electrons and 
Nd down electrons, we need to specify positions of up electrons as Ua,m ("^ = 1, ■ ' ' ^ ^u) 
and that of down electrons as da^n = l,''';^^)- They satisfy 1 < Ua,i < Ua,2 < 
■ ■ ■ < Ua,Nu < ^ 1 < da,i < da,2 < ■ ■ ■ < da,Na < ^- the present case, we have a 
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half-filled Hubbard model whose ground state is given with = = N/2. 

Note that for a pair of different k points k ^ k\ (^'^\c\ ^cy = {k,(y\k',a) = 0. 
The charge density is thus represented as, 

n(r)= ^(vI/|V;t(r)V^,(r)|v]>) 

cr k,k' 
cr k 

= J2J2\Mr)\Mk,a). (21) 

a k 

n{k,a) is the momentum distribution given by |\E'). The system is found in a 
paramagnetic state and 0fc(r) and n{k, a) = n{p) lose the spin dependence, in which the 
crystal momentum p = 2Tik/ (Na) is used. 

In this simulation, the value of U is approximated to be f/ = 5.2to, which is roughly 
estimated by the result of the Hydrogen molecule in the last sub-section. In the obtained 
self-consistent solution, the transfer terms tij are given by the Fourier transformation 
of the Kohn-Sham eigenvalues e{n,p). Only the Is band (n = 1) is used to construct 
tij ■ 

We show the single-particle dispersion of Eq. (ITT!) and the momentum distribution 
of the obtained self-consistent solution in Figures 4 and 5, respectively. The many-body 
model Eq. f|T2l) becomes a kind of the one-dimensional Hubbard model. We can see 
that ?T,(r) is almost unchanged by introduction of U, while (n^) is suppressed by the U 
term. This is seen in the dispersion relation of e{n,p), which is almost the same for cases 
with a finite U and without U. On the other hand, when U = 5.2to, n{p) is completely 
different from that of the free Fermion. The feature of n{p) as well as the filling factor 
of the system tells that the system is in a Mott insulating phase. 

7. Discussion 

We have a concept of the fixed-point Hamiltonian in our theory, which is defined in the 
whole phase space of the original problem. This fact is in contrast to the usual idea of the 
renormalization group. The smearing process in our formulation is the self-consistency 
loop, in which effective interaction processes are rebuilt via the redefinition of 0i(r). 
On the contrary to the usual renormalization group analysis, in which the zooming out 
process inevitably smearing out microscopic details of the order parameter, the central 
order parameter n(r) is kept its microscopic structure in the present formulation of MR- 
DFT. A reason why we can reconstruct the effective many-body Hamiltonian comes from 
the flexibility of EKSS based on the density functional theory. 

In the present formulation of EKSS, people might think that the reference 
calculation is inevitable to obtain the value of U. If we utilize LDA for the residual 
exchange- correlation energy functional, the approach may seem close to established 
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Figure 4. The Kohn-Sham eigen values e{n,p) of Eq. [11] which gives the single- 
particle dispersion of a Hydrogen chain. The value ofU/to — (crosses) or 5.2^0 with 
to (pluses) being the transfer energy between neighboring atoms. 




Figure 5. The Fermion momentum distribution n{p) of the Hydrogen chain with 
iV = 10 atoms. The value of U/Iq — (crosses) or 5.2io with to (pluses) being the 
transfer energy between neighboring atoms. 
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Figure 6. The total energy i5tot for the Hydrogen molecule with the inter-atomic 
distance R[A] obtained by EKSS. Depending on the fluctuation (nf), -Etot increases 
monotonically. 



LDA+U. Now, we will propose an indicator to find out the clue of change in 
the fluctuation appearing in the system. We also discuss a method to detect the 
Mott insulating transition in MR-DFT. Due to these characteristic factors, EKSS is 
qualitatively and quantitatively different from the known LDA+U approaches. 



7.1. An indicator for fluctuation suppression 

We analyze the EKSS result of the Hydrogen molecule to test the refinement of the 
residual exchange- correlation energy functional. In Figures 6, 7 and 8, we show a total 
energy, the kinetic energy and the Hartree term of the system. Here, the definition of 
the total energy is. 



Etot = (^|T|*) + 



1 f r2^(r)r?,^(r' 



■d^rd^r 



3„' 



(22) 



in which contribution of the U term is omitted. |\E') is obtained by solving Eq. (fT2l) . so 
that the state is a correlated Fermion state. The kinetic energy and the Hartree term 
denote Ekin = (^|T|^) and 



E- 



Hartree 



■d-'rd^r' 



Now, we have another expression for E^. Consider the minimizing I^E') of 
(^'|T + Kedl^') which gives ncslr) and is the solution of Eq. f|T2|) . Then, we have. 



■d'^rd^r' 
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Figure 7. The kinetic energy E'kin for the Hydrogen molecule with the inter-atomic 
distance R[A] obtained by EKSS. The value is written as a function of which is 
controlled by U. 
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Figure 8. The Hartree term i^Hartrcc for the Hydrogen molecule with the inter-atomic 
distance R[A] obtained by EKSS. The value is written as a function of which is 
controlled by U. 
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+ Er^n^] + j d^rvcxtir)n^{r) 
= m + V^edlvP) + I J ^^^^^^d^d'r' + E^^M 

+ mill ($'|f 1$') - mill (^'|f + Kcdl^') + / c/Vfext(r)n^(r) 

Z J |r — r I 

+ / d^rve^t{r)n^{r) + min ($'|r|<l>') - (^|f |^). (23) 

J $'^n>j 

Thus we may write Eo as, 

Eo = Etot + mill ($'|f 1$') - (^|r|$). (24) 

This is another exact expression of the true total energy of the electron system. Note 
that the U term does not appear in the formula, although it affects |\E'). People might 
find that the above expression can be used to avoid the double counting problem. Let us 
evaluate Eq within the approximation utilized in Sec. [61 Now, look at the kinetic energy 
i?kin for the Hydrogen molecule. (Figure 7) When R < 1.0 [A], the value decreases with 
decreasing (n^), which is controlled by increasing U. Namely, the horizontal axis is the 
direction of increasing U. This reduction in the kinetic energy is caused by expansion 
of the wavefunction in the real space. Actually, the Hartree term decreases and the 
electron-ion potential terms reduces their absolute values. In this range, as seen in the 
shift in the Hartree term, n(r) expands with increasing U. Thus, we find that i^^km 
decreases in a weakly correlated regime {R < 1.0 [A]) by increasing U. 

Let us compare the result with the cases with R > 1.0 [A]. In this region, i?kin 
increases by increasing U. If we look at -^Hartree, we see that the value does not change 
so much and is almost constant, when R > 2.0 [A]. This fact means that n(r) is nearly 
unchanged. What the U term does in this regime is that it only shift the internal 
fluctuation. Thus, the value of i^kin increases. Now look at the expression of Eq. (^^. 
The true value of Eq is estimated by adding the kinetic energy of an uncorrelated 
Fermion system ($'|T|$'), which has n{r) = n^(r), and subtracting (\E'|T|\E^) from -Etot- 
In this example, since n(r) is nearly unchanged against shift in U for R ^ 1.0 [A], 
min$/^„^($'|T|$') for finite U may be approximated by E'km ioT U = 0. The result 
suggests that Eq is almost unchanged by increasing U, while the state l^') becomes a 
correlated state. 

On the other hand, if we detect decrease of -Ekin by increasing U, this suggests 
that minimizing $' of ($'|T|$') should be close to \E'. We have an inequality, 
min$/^„^($'|T|$') < (\1'|T|\E'). Thus, Eq evaluated for finite U is nearly the same as 
i^tot- However, i^tot increases by introduction of U. When we have the weakly correlated 
regime R < 1.0 [A], the U term is not necessary for the proper description of the system. 

As a result, we conclude that we can utilize [/-dependence of i^km = (^E^ITI^^) to 
detect occurrence of the Coulomb suppression in a correlated electron system. Once we 



A first-principles calculation for correlated electron systems 



21 



have a properly designed method to estimate min$/^„^ ($'|T|$'), EKSS works as a first- 
principles calculation method for the correlated electron system in general even without 
a reference calculation prepared for each individual problem. The target systems for 
EKSS include the Mott insulating state. Actually, we know a numerical algorithm [38] 
to obtain the Legendre transform, 

E{n) = sup 

V 

7.2. A test for the Mott insulator 

To test the conduction property of the system within DFT, we may be able to utilize the 
next technique of the momentum boost. Let us consider a twisted boundary condition 
for our simulation. 

^(r + L^e^) =exp(i0)^(r), ^(r + L^e^^) = ^(r), ^(r + L,e,) = ^'(r), 

Bj and Li {i = x, y, z) are the unit vectors and the length of a simulation cell. 
The density- functional theory holds for any fixed 9. Let us shift 9 from zero to 2tt 
adiabatically and obtain the lowest energy eigen value Eq{9). Then we can connect 
Eo{9) and draw a graph of Eo{9) function of 9. 

According to the Kohn argument, [39] we can identify the Mott insulating state 
by looking at the period of Eo{9), although we may see only the lowest edge of the 
whole Eq{9). If formation of a gap in the flow of Eq{9) is detected by changing the 
lattice constant, for example, the system undergoes the Mott transition. Actually, a 
complete test using the one-dimensional Hubbard model showed the period 27r for the 
half-filled band, that is useful for the characterization of the ground state. [IQ] If the 
system is described in the Kohn-Sham scheme with LDA, however, the period would 
not change from the value of a metallic state. This failure would be recovered by the 
introduction of the U term in the Kohn-Sham scheme. If we ask the system to reproduce 
the local fluctuation, modification of the Kohn-Sham system naturally makes the system 
interacting. This is a way to model the stiffness of the Mott insulating state against 

Q 

the boost induced by the imaginative magnetic flux, which amounts to — $o with the 

27r 

unit flux $0- The nature of the ground state is modified via a change in the charge 
fluctuation. 

Acknowledgement 

The authors are grateful for many occasions of discussion with many researchers. K.K. 
especially thanks professors, H. Aoki, S. Tsuneyuki, H. Kamimura, M. Tsukada, M. 
Imada, M. Ogata, A. Hasegawa, H. Akai, H. Katayama-Yoshida, H. Kasai, M. Higuchi, 
K. Higuchi, Y. Morikawa, J. Yamauchi, and doctors, T. Ogitsu, K. Kobayashi and, Mr. 
M. Takahashi. This work was supported by a Grant-in-Aid for Scientific Research in 
Priority Areas "Development of New Quantum Simulators and Quantum Design" (No. 
17064006), the 21st century COE program "Core Research and Advanced Eduation 



min(\E'| 



T + / drv{r) (n(r) - n{r)) > |\E') 



A first-principles calculation for correlated electron systems 



22 



Center for Materials Science and Nano Engineering", Grants-in-Aid for Scientific 
Research (No. 15GS0213) and also by a Computational Nanoscience program "Grid 
Application Research in Nanoscience-National Research Grid Initiative (NAREGI)" of 
the Ministry of Education, Culture, Sports, Science, and Technology, Japan. 

References 

[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 

[2] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 

[3] M. Levy, Proc. Natl. Acad. Sci. (U.S.A.) 76, 6062 (1979). 

[4] M. Levy, Phys. Rev. A 26, 1200 (1982). 

[5] E. Liob, Int. .J. Quantum. Chcm. 24, 243 (1983). 

[6] N. Hadjisavvas and A. Theophilou, Phys. Rev. A 30, 2183 (1984). 

[7] E.S. Kryachko and E.V. Ludena, Energy Density Functional Theory of Many-Electron Systems, 
Kluwer, (1986). 

[8] R.G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Plenum, (1989). 

[9] R.M. Dreizler and K.U. Gross, Density Functional Theory, Springer, (1990). 
[10] L.F. Mattheiss, Phys. Rev. Lett. 58, 1028 (1987). 
[11] J.Y. Freeman and J.-H. Xu, Phys. Rev. Lett. 58, 1035 (1987). 

[12] K. Takagahara and H. Harima and A. Yanase, Jpn. J. Appl. Phys. 26, L352 (1987). 
[13] T. Oguchi, Jpn. J. Appl. Phys. 26, L417 (1987). 

[14] W.E. Pickett and H. Krakauer and D.A. Papaconstantopoulos and L.L. Boyer, Phys. Rev. B. 35, 
7252 (1987). 

[15] R.V. Kasowski and W.Y. Hsu and F. Herman, Solid State Commun. 63, 1077 (1987). 
[16] R.V. Kasowski and W.Y. Hsu and F. Herman, Phys. Rev. B 36, 7248 (1987). 
[17] J. Hubbard Proc. Roy. Soc. (London) A276, 238 (1963). 

[18] D. Baeriswyl and D.K. Campbell and J.M.P. Carmelo and F. Guinea and E. Louis, The Hubabrd 

Model, Plenum, (1995). 
[19] K. Yamaguchi, Chem. Phys. Lett. 68, 477 (1979). 
[20] G.C. Lie and E. Clementi, .J. Chem. Phys. 60, 1275, 1288 (1974). 
[21] F. Moscardo and E. San-Fabian, Phys. Rev. A 44, 1549 (1991). 
[22] B. Miehlich and H. Stoll and A. Savin, Mol. Phys. 91, 527 (1997). 
[23] J. Grafenstein and D. Cremer, Chem. Phys. Lett. 316, 316 (2000). 
[24] J. Grafenstein and D. Cremer, Mol. Phys. 103, 279 (2005). 

[25] S. Yamanaka and K. Nakata and T. Takada and K. Kusakabe and J.M. Ugalde and K. Yamaguchi, 

Chem. Lett. 35, 242 (2006). 
[26] K. Kusakabe, J. Phys. Soc. Jpn 70, 2038 (2001). 
[27] K. Kusakabe, cond-mat/0505703 (2005). 

[28] V.L Anizimov and J. Zaanen and O.K. Andersen, Phys. Rev. B 44, 943 (1991). 
[29] I.V. Solovyev and P.H. Dederichs and V.I. Anisimov, Phys. Rev. B 50, 16861 (1994). 
[30] A.I. Liechtenstein and V.I. Anisimov and J. Zaanen, Phys. Rev. B 52, R5467 (1995). 
[31] J.E. Harriman, Phys. Rev. A 6, 680 (1981). 
[32] G.H. Wannier, Phys. Rev. 52, 191 (1937). 
[33] N. Marzari and D. Vanderbih, Phys. Rev. B 56, 12847 (1997). 
[34] K. Kusakabe and M. Takahashi and N. Suzuki, Physica B, 378-380, 271 (2006). 
[35] M. Takahashi and K. Kusakabe and N. Suzuki, Proc. 28th Int. Conf. PHys. Semicon. (Wien), in 
press. 

[36] J.F. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 
[37] N. TrouUier and J.L. Martins, Phys. Rev. B 43, 1993 (1991). 
[38] K. Kusakabe, unpublished. 



A first-principles calculation for correlated electron syst 

[39] W. Kohn, Phys. Rev. 133, A171 (1964). 

[40] K. Kusakabe, J. Phys. Soc. Jpn. 66, 2075 (1997). 



